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Abstract 

Real data are constrained to finite sampling rates, which calls for a suitable mathematical de- 
scription of the corrections to the finite-time estimations of the dynamic equations. Often in the 
literature, lower order discrete time approximations of the modeling diffusion processes are con- 
sidered. On the other hand, there is a lack of simple estimating procedures based on higher order 
approximations. For standard diffusion models, that include additive and multiplicative noise com- 
ponents, we obtain the exact corrections to the empirical finite-time drift and diffusion coefficients, 
based on Ito- Taylor expansions. These results allow to reconstruct the real hidden coefficients from 
the empirical estimates. We also derive higher-order finite-time expressions for the third and fourth 
conditional moments, that furnish extra theoretical checks for that class of diffusive models. The 
theoretical predictions are compared with the numerical outcomes of some representative artificial 
time-series. 
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I. INTRODUCTION 



Many fluctuating random phenomena, including turbulent diffusion, polymer dynamics 
or asset price evolution, can be modeled by an univariate Ito-stochastic differential equation 
(SDE) of the form bellow, characterizing a diffusive model: 

dX t = D 1 (X t )dt + y/2D 2 (X t )dW t , (1) 

where W t is a Wiener process, Di(X t ) is the coefficient of the slowly varying component 
(called drift coefficient) and D 2 (X t ) is the coefficient of the rapid one (called diffusion coef- 
ficient). 

For sufficiently smooth and bounded drift and diffusion coefficients, the associated prob- 
ability density function (PDF) P(x, t) = P(X t = x, t) is governed by the corresponding 
Fokker-Planck equation l| 

d t P(x,t) = -d x [D 1 (x)P(x,t)]+d xx [D 2 (x)P(x,t)]. (2) 

Here, we are concerned with the empirical access to unknown drift and diffusion coeffi- 
cients of stochastic processes. For an ideal time series X t generated by Eq. ([I]) and sampled 
with a sufficiently high resolution on a long time period, the original coefficients can be 
perfectly reconstructed. For stationary processes, the coefficients D k (x), with k = 1,2 can 
be directly estimated from the conditional moments l] as: 

D k {x) = ]xmD k (x,r) , (3) 

T— >0 

where 

D k {x,r) = ^([X t+ r - X t f)\ Xt=x , (4) 

with (•••) denoting statistical average and \x t =x meaning that at time t the stochastic 
variable assumes the value x. 

Conversely, for general Markovian stochastic processes, the time evolution of PDFs is 
governed by a generalization of Eq. (121), namely 

d t P{x, t) = Y,(-d x ) k [D k {x)P{x, t)\ . (5) 

fc>0 

with coefficients D k (x) given by Eqs. ©-(jl]), for any integer k > 1. For diffusive processes, 
Eq. (jSJ) reduces to Eq. (j2J). Therefore, processes governed by the Ito-Langevin Eq. (TTJ) must 



11] 



furnish null coefficients Dk, for k > 3. Pawula theorem |lj simplifies this task by stating 
that if -D4 is null, all other coefficients with k > 3 are null as well. The coefficient -D4 is then 
a key coefficient to be investigated, in order to establish the validity of the modeling of data 
series by Eq. 

However, due to the finite sampling rate of real data, numerical estimations of can 
not always be straightforwardly extrapolated to the limit r — > in Eq. ([3]). In such cases, 
one accesses only the finite-r estimation of the coefficients given by Eq. (j4]), which may 
significantly differ from the true coefficients D^. This is specially relevant when r is large 
compared to the characteristic timescales of the process. 

Some authors 3] have introduced finite sampling rate corrections to the coefficients 
Di(x, r) and D 2 (x, r), by deriving expansions for the conditional moments up to some speci- 
fied low order of r, directly from the Fokker-Planck equations. Applications of this approach 
have already been implemented for those coefficients up to second order j^]. The error in the 
finite-r estimated coefficients Dk can also be derived from the stochastic Ito- Taylor expan- 
sion [4] of the integrated form of Eq. (JTJ). Within this line, the first order expansion of drift 



and diffusion coefficients was recently presented in Ref. [5(. However, low order corrections 
may be inappropriate when the convergence of the limit in Eq. ([3]) is slow j^, 6]. Moreover, 
there is no a priori knowledge of whether the sampling rate is fine enough to justify the use 
of the lowest order approximation. 

In the present work, we investigate those issues for diffusion models defined by Eq. (JTJ 
with linear drift coefficient, namely, -Di(x) = — aix, representing an harmonic restoring 
mechanism, and quadratic state-dependent diffusion coefficient, namely, D 2 (x) = b + b 2 x 2 . 
This class encompasses some of the most common models of the theoretical literature. In 
fact, this equation is frequently found in a diversity of processes, from turbulence to fi- 
nance M, M . Moreover, the obtained results are also valid for another class of SDEs with 
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additive-multiplicative noises [9|, H0[ , given by 

dX t = -axXtdt + V%>odWt + ^h 2 X t dW[ , (6) 

where W t , W[ are uncorrelated Wiener processes. 

For discretely sampled data at intervals r, we will derive, from the stochastic Ito- Taylor 
expansion, finite-r expressions for the parameters {01,60,^2}, up to infinite order. These 
exact expressions will allow us to reconstruct the true drift and diffusion coefficients from 



their empirical finite-time estimates. As a corollary, one can determine up to which value of 
t a given order of truncation is reliable (within a fixed tolerance), or reciprocally, which is 
the sufficient order for a given r. 

Furthermore, as empirical estimates suffer from finite-r effects, one always gets non- 
null -D4. Therefore, the evaluation of the corrections for this coefficient is crucial for a 
suitable probe of the diffusive modeling. In this work, we also derive finite-r expressions for 
coefficients D3 and D4, which furnish extra theoretical tests of consistency for the diffusive 
models considered. 

Our theoretical findings are corroborated by the outcomes of exemplary artificial time- 
series generated by Eq. ([TJ. 

II. EXACT CORRECTIONS FOR DRIFT AND DIFFUSION COEFFICIENTS 

Let us consider the Ito formula 4j, for a given function F of the stochastic variable X t 

dF = (d t F + D^xF + D 2 d xx F)dt + ^2D 2 d x FdW 

= L Q Fdb + l}FdW, (7) 

and its integrated form 

/t+T pt+T 
L°F(X s )ds + J L l F(X s )dW s . (8) 

Let r be the sampling interval of state space observations. By applying Ito formula 



to the functions Di(X s ) and y 2D2(X S ) in the integral form of Eq. (CQ): 

X t+T = X t + [ + D 1 {X s )ds+ f + y / 2D 2 {X s )dW s , (9) 
Jt Jt 



one finds 



X t+T = X t + £ + (d^XJ+J L°D 1 (X' s )ds' + J L 1 D 1 {X sl )dW s ^jds 

/ 2D 2 (X t ) + jT L°y/2D 2 (X' s )ds' + jT LW2D 2 (X s ,)dW s AdW s . (10) 



After iterated applications of Ito formula, one gets an expression for the increment of the 
stochastic variable in terms of multiple stochastic integrals 4j: 

X t+T -X t = J2ca k (Di,D 2 )I ak , (11) 



where a k = {juh, ■ ■ ■ ,3k), with ji = 0,1 for all i, c ak (D 1 ,D 2 ) = 
LP 1 LP 2 . . . LP k ~ 1 LP k and I ak are multiple stochastic integrals of the form I ak = 

It +T It t+tk !t +tk ' X ■ ■ ■ It t+t2 dt i ■ ■ ■ K-IK, with dt i = dU and dtj = dW t . 

By inserting Eq. ffTT]) into Eq. (j3J) and performing the averaging, for k = 1, 2, we achieve 
analytical expressions for the fmite-r drift and diffusion coefficients, up to arbitrary order 
in powers of r. The resulting expressions preserve the linear and quadratic x-dependence, 
respectively and can be written as: 

Di(x, t) = —ai{r)x 

D 2 (x,r) = ~b (r) +b 2 (r)x 2 . (12) 

Hence, we are led to the theoretical relation between the finite-r coefficients {ai, 60,62} 
and the true ones {ai, 60, b 2 }, namely, 



j>0 u ^ ; ' 



j>0 U ^ ; ' 

62 ( r ) = E ^riv • ( } 

Details of the derivation of Eqs. (!TH]) - (fT5|) can be found in the Appendix. 

By restricting the expansions (1X31) - (TT5T) to some common finite power n, one gets the 
respective nth-order approximation. This result extends previous findings of first [5] and 
second jjj order terms. 

Notice that Eq. (TX5T) is uncoupled, meaning that the estimated harmonic stiffness a\ 
is not affected by the exact noise components. Moreover, from Eq. (115j) . the estimated 
multiplicative noise parameter b 2 does not depend on the exact additive noise component. 

Summing the series in Eqs. ffT3l) - ffT5l) up to infinite order, and defining Z = exp(— air) 
and W = exp(— 2b 2 r), we find the exact finite-r expressions: 

1 - Z 
a\ = , 

T 

~ b 1 - Z 2 W 

b ° = 1 n ' 

ai — b 2 It 

~ l-Z 1-Z 2 W 

b 2 = . (16) 

t It 



Notice that lim r _».o{Si, bo, 62} = {^i, 60,62} holds. 

From Eqs. (|T6l) . we obtain an invariant relation among the estimated and exact parame- 
ters, namely, 

Qi - b~2 _ a\ - b 2 ^ 
b b 

The meaning of this invariance can be drawn, for instance, from the stationary PDF 
P*(x) associated to the corresponding Fokker-Planck equation given by Eq. ([2]). With the 
present choice of drift and diffusion coefficients, for ai, 60 > 0, b 2 > 0, one has: 

P*(x) = P /[l + % 2 }^ + \ (18) 

with P Q a normalization constant. This solution is of the g-Gaussian form Q], for which, if 
a i ~ b 2 > 0, the variance is finite with value a 2 = &o/( a i — b 2 ). Hence, Eq. (flTl) represents the 
uphold of the data variance under changes of sampling intervals. For b 2 = 0, one recovers the 
Gaussian stationary solution and its variance relation. Notice also that, from Eqs. (|T3l) - (|T5l) . 
Eq. ( flTl) still holds if one considers partial corrections of the parameters up to any common 
order n of truncation of the sums. 

Let us remark that the results presented in Eqs. (|T2l)-(|T6l) are valid even when the variance 
is infinite. However, we will deal only with finite variance cases and consider normalized 
data (with unitary variance), which only implies a rescaling of b — > bo /a 2 . Then, 

ai = b + b 2 . (19) 

Taking into account the constraint ( fl9l) . from Eqs. (fl6l) . the exact finite-r expressions 
are: 

a, = 1 " e *p(-" T ) ( 2o) 

60 = 1 " CTP 2" 26 ° T) . (21) 
Eqs. (1201) and (12T1) can be readily inverted to extract the true parameters from their 



finite-r estimates: 

ai = in(i : ° it)) , (22) 

bo = In(l-W)) p3) 
Notice that a^r (and also 2b r) can not be greater than unit. 



In what follows, we fix the timescale r = 1. A different choice would simply lead to a 
rescaling of the parameters (ai, 60, b 2 ) — > (rai, rb , r& 2 ) ■ 




FIG. 1: Drift and diffusion coefficients for the O-U process. Symbols correspond to the numerical 
computation for artificial series (10 5 data), synthetized with the values of a\ = bo (62 = 0, in accord 
with constraint indicated on each panel. Lines represent the coefficients given by Eqs. (|12|) . 

using the theoretical r-expansions fll3])-(|15|). at different orders of truncation. The darker the color, 
the higher the order, from first up to fifth order. The infinite order (exact expression) is represented 
in thick black lines. The zeroth order, corresponding to the true values, is plotted in dashed lines. 



Now we investigate the importance of finite-r effects for discretely sampled realizations 
of representative known diffusive processes. To this end, we generated artificial time-series 
through numerical integration of Eq. (pQ), by means of an Euler algorithm with timestep 
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dt = 10~ 3 , recording the data at each 1/dt timesteps, in accord with our choice r = 1. Our 
theoretical results for D\ and D 2 will be compared to the ones numerically computed from 
the time-series, through Eq. (j3J). 

The particular case b 2 = 0, corresponding to the Orstein-Uhlenbeck (O-U) process and 
the general case with multiplicative component 62 > will be investigated separately. Fig. [T] 
shows the results for the artificial series with known values of the parameters a\ = bo (b 2 = 0), 
together with our theoretical predictions. The exact theoretical expressions reproduce the 
numerical (finite-time) outcomes. Comparing the panels in Fig. [TJ it is clear that, the 
larger a%, the slower the convergence to the observed coefficients. The results for a\ > 1 also 
illustrate the entanglement one may find in large-r measurements, specifically, an oscillatory 
convergence of d\ and an alternating signal of b 2 . 

In Fig. [21 we plot the numerical computations for artificial time-series together with 
analytical predictions for b 2 > 0. Again, the theoretical approximations present slower 
convergence as a\ increases while the exact theoretical expressions agree with finite-time 
estimates directly obtained from the time-series. Moreover, the actual value of b 2 sets the 
convergence rate of b 2 . 

All these results raise a question about the domain of validity of lower order approxima- 
tions presented before in the literature. Let us investigate this issue quantitatively. Given 
di, obtained from numerical (finite-time) evaluation, the exact value of a\ can be recovered 
from Eq. (1221) . Approximate values 0^ can be obtained by inversion of Eq. ( fT3l) truncated 
at order n. Figure ([3]) illustrates function of n, for different values of d\. Clearly, 

convergence to the true value a\ is attained (within a given tolerance), at different orders 
that depend on the value of d\. For instance, for d\ > 0.5, an order larger than two is 
required. Convergence is faster for smaller d\, that is, as soon as \jd\ becomes large com- 
pared to the timescale r = 1. For bn we obtained a very similar convergence scheme (not 

nn 

shown) . In Refs. [3|, |5[ , the fitness of low order expressions for O-U processes results from 
the particular employment of d\T < 0.5. However, this may not be the case when dealing 
with generic empirical data. 
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X 



FIG. 2: Drift and diffusion coefficients for the general process with multiplicative noise. Symbols 
correspond to the numerical computation for artificial series (10 6 data), synthetized with the values 
of ai (= 60 + 62 1 m accord with constraint (|19p ) and 62 indicated on each panel. Lines are as in 

Fig.m 
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FIG. 3: Dependence of on the order n of the approximation given by Eq. (|13p . for different 
values of a\ (panel (a)). Dotted lines correspond to the respective true values {a{) and missing 
points denote the absence of real solutions. Panel (b) exhibits the order at which the limiting value 
is attained (within 5%) as a function of a\. 

III. HIGHER-ORDER COEFFICIENTS 



Inserting the Ito- Taylor expansion Eq. (Hip into Eq. (J3J) and performing the average for 
k — 3, 4, we also computed the fmite-r expansion for D 3 (x,t) and D 4 (x,r). The resulting 
expressions are invariant functions of x, namely: 



D 3 (x,t) = -ci(t) + c 3 (t)x 3 
Z) 4 (x,t) = d (r) + d 2 (T)x 2 + (i 4 (r)x 4 



(24) 
(25) 



For the particular case h 2 = 0, we were able to derive the infinite order expansion for the 
r-parameters: 
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Ci(T 

J (r 
d 4 (r 
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E- 

i>o 
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+1 



1 



2 (j + 1)! 
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(26) 



4; _ 2 J [_ 0l] i-i 



i>o 



j>o 
4J _ 3^+1 



2 (j + 1)! ' 
2 x 4 J - 3 J+1 + 1 [-ai] j 



2 + 



6 



(j + i)! 



(27) 



Notice that, lim T _, {ci, c 3 , d , o?2, ^4} = , as expected, and that the relevant parameter for 
the rate of series convergence is a\. Summing the series (J26l) -( l271) up to infinite order, and 
recalling that Z = exp(— g^t), one obtains 

b (l-Z)\l + Z) 



d 2 (r 



ai 2r 




{1-Zf 




6t ' 




hi (i-znu 


-z? 


a\ 8t 


■ 


b (1 - Zf{l 4 


-z) 



(28) 



di 4r 

(1 - zf 

24r 



(29) 



Fig. H] shows the numerical computation of D 3 (x,t) and D±(x,t) for the same artificial 
series as in Fig. [TJ For comparison, the theoretical estimates at different orders of truncation 
of the series in Eqs. (j26"l) - (l2"7l) are shown. Notice that, although D 3 , = for the diffusive 
processes considered here, their finite-time counterparts have cubic and quadratic forms. 
Indeed, the exact theoretical expressions given by Eqs. (j28]) -(129~ ]) reproduce the numerical 
outcomes, validating our approach as furnishing meaningful tests for O-U models. However, 
for di > 1, rich pictures for the low order approximations of D 3 and D4 arise, which hinder 
the asymptotic estimation. 

For the general case with 6 2 > 0, we computed the third-order r-expansions for D 3 (x, r) 



and D^{x,t). Each power of r of order j < 3 has pre-factor denoted by £)g and D 
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FIG. 4: Third and fourth order coefficients for the O-U process. Symbols correspond to the 
numerical computation for the same artificial series of Fig. [TJ Lines represent the coefficients given 
by Eqs. (|24p - (|25p . using the theoretical r-expansions (|26p - (|27p . at different orders of truncation. 
Colors as in Fig. [TJ 

respectively. We find, for Z) 3 : 
D® = 0, 

L>3 = — b ax — b 2 ax , 

Z)f = -b a(9ai — 16b 2 )x a(a\ — 13a!& 2 + 166n)^ 3 , 

6 6 

= -i& a 2 (8a! - 13b 2 )x + — a(3a? - 32a\b 2 + 7A ai b 2 2 - 52b 3 2 )x 3 , (30) 
6 12 

where a = a\ — 26 2 . At all orders, _D 3 vanishes if ai = 2b 2 . 



12 




FIG. 5: Third and fourth order coefficients for the general process with multiplicative noise. Sym- 
bols correspond to the numerical computation for the same artificial series of Fig. [21 Lines represent 
the coefficients given by Eqs. (|24p -(|25 p . using the theoretical r-expansions ([30]) - ([3Tj) . at different 
orders of truncation up to third order. Colors as in previous figures. 
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For D4, we have: 

Df ] = 0, 

= \bl + b b 2 x 2 + \f 2 x' , 

Df = -77&o(3ai - 7b 2 ) + h (3al - 30 ai b 2 + 52b 2 2 )x 2 + \b 2 {?>al - 2A ai b 2 + 38bj)x 4 , 
6 b o 

I)f = -blClal - 34a x 6 2 + 436?) - -b (6al - 65a?6 2 + 206^^ - 206b 3 2 )x 2 
6 6 

+ - 36a ? & 2 + 21ba\b\ - 736ai^ + 652^)x 4 . (31) 

In Fig. El we show the numerical computation of D3 and -D4 for the same artificial series 
as in Fig. According to our theoretical results, a\ = 2b 2 is a threshold between positive 
and negative slopes of D 3 , as illustrated in Fig. El For comparison, we also show up to the 
third order theoretical estimates Ylj>i ^3 an< ^ 2~2j>i according to Eqs. fl30]) - fl3T|) . 

Coefficientes D 3 and D4 provide further tests of validity of the diffusive modeling. From 
Fig. El third order estimates furnish suitable forecast of the numerical finite-r measurements 
for small enough a±. In such cases, once obtained ai, 60 from Eqs. (I22l) -( l23l) . these values 
can be used in theoretical equations for D3 and -D4, to check if the corresponding non-null 
coefficients can be attributed to finite-r effects. 



IV. SUMMARY AND FINAL COMMENTS 



For an important class of diffusion models with additive-multiplicative noise, we have 
derived exact formulas that connect the empirical discrete-time estimates with the actual 
values of the parameters of drift and diffusion coefficients. Additionally, we also provided 
theoretical expressions for higher-order coefficients which serve as a further probe for the 
validity of this class of diffusive models. 

Our results allow to access the generating stochastic process. A possible procedure to 
identify it and its parameters can be summarized as follows. When numerical computation 
of the coefficients from a real timeseries yields linear and quadratic forms for Di and D 2 
(which is a frequent outcome), the values of {01,60^2} can be obtained from curve fitting. 
The present model (with or without multiplicative component) would be adequate when i) 
di < 1 and ii) (ai ~ 6 + b 2 ). For d\ < 0.5 the second order correction would be enough 
to recover aj, otherwise larger order corrections should be considered. Once d\ and b are 
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known, Eqs. (|22|) and (1231) . allow to obtain, exactly, the original parameters ai,bo, hence 
also b 2 = a\ — b . Values such that b 2 << «i,&o (hence, a± ~ 6 ) point to a simple O-U 
process, otherwise a multiplicative term may be also present. In both further check 

consists in the analysis of higher order coefficients, e.g., to see whether a non-null -D4 can 
be attributed to finite-time corrections. 

By analyzing the O-U process, we also found that a low sampling rate would significantly 
affect the diffusion coefficient estimate, by adding an extra quadratic term. Thus, the 
detection of a quadratic D 2 does not imply the existence of multiplicative components in 
the actual process. Moreover, estimations of D3 and D4 from low-order approximations 
would lead to results inconsistent with the empirical outcomes. 

The obtained formulas also allow to quantify the errors induced by a finite sampling 
rate r in the numerically estimated coefficients. The analytical results indicate that, in 
order to grasp the true values of the parameters from the knowledge of the observed ones, 
the required correction depends strongly on the (hidden) inverse time a±. Our work shows 
that one should be careful when applying low-order finite-r corrections for diffusion models. 
Furthermore, as shown in Fig. [3j our results provide a criterion, from the knowledge of ai, 
to determine the required order n, or equivalently, up to which value of r the respective 
approximation is reliable. 
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Appendix 

Rewriting Eq. (J4]) according to the notation introduced in Eq. (TTTj) . one has 

TD 1 = {AX) = Y,c*>(D l ,D 2 )(I ak ). 

&k 

Only multiple stochastic integrals I ak such that at = (0, 0)^ have non-null average, being 
(I(o,...,o) k ) = r k /k\. From the iterated application of ltd formula to Eq. (jHJ), C( 0j ... i o) fc (-Di, D 2 ) = 
(L°) fc_1 Di. These are general results independent of the particular form of D\ and D 2 . 
Noticing that, from Eq. ([7]), L° = d t + Did x + D 2 d xx , and that Di is time-independent and 
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linear in x, then, c Qfe = Di(D[) k 1 = (— a\) k x. Finally, 

D 1 = (AX}/T = J2y(-*i) k r k - 1 x, 

k>l 

which is of the same functional form of the true Di and can be identified with —aix, so that 

*!>i 

which gives Eq. ffTB"]) . 

For the second conditional moment, one has 

2rD 2 = ((AX) 2 ) = ^ c an c Pm (I a Jp m ). (32) 

From the definition of c ak in Eq. (fTTj) . if 62 = (then D2 is constant), only two classes of 
terms in Eq. (1321) are non-null, those with: 

i) a n = (0, 0) n and (3 m = (0, 0) m and 

ii) a n = (l,0,...,0) n and (3 m = (1, 0, 0) m . In those cases, the products c an cp m take the 
values: 

i) {D l f{D[) k - 1 = (-ai) fc+1 x 2 (with k = m + n-l), 

ii) 2D 2 (D[) k = 2b (-a 1 ) k (with k = m + n - 2). 

In order to evaluate the averages of products of multiple stochastic integrals, it is useful 

to recall that (/(o,... 1 o)„J(o,...,o)J = and that ( J (i,..,o) J(i,...,o)J = {n+m -\)(n-iy.( m -iy. El- 
Then, summing over all the pairs (n, m) contributing to the order r fc+1 , one obtains: 

i) Y(I U )/r k+1 = 1 Y k ( k+1 ) = 2 2k ~ 1 

ii) J2(Ia n If3 m ) /r k+1 = ( fc+1 )| Z^=0 (j) = (fc+1)! ■ 

Finally, from Eq. (I32l) . we arrive at 

D 2 = <(AX) 2 )/(2r) = \ g ^^^(-oO^V + (^^o(-«i) fc ) T fc , (33) 

which can be cast in the form b 2 x 2 + allowing to identify b 2 and b Q with functions of the 
true parameters, as 

k>0 K ' 

^ = E(^T(-«0 t+ V'. (34) 

fc>0 v ; ' 
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For the general case 62 > 0, a similar but tricky derivation leads to Eqs. f|T4l) - (|T5l) that 
generalize the expressions (1341 . 

Proceeding with the third order, products of three multiple integrals appear. For 6 2 = 0, 
there are two types of products I oln Ip m I ll contributing to ((AX) 3 ): 

i) a n = (0, 0) n , p m = (0, 0) m , 7; = (0, 0)j and 

ii) a n = (0, ...,0) n , p m = (1, 0, 0) m , 7; = (1, 0, 0)/, with pre-factors proportional to 
(£>i) 3 (£>;) fe - 2 = (-ai) fc+1 x 3 and to 2D 2 D 1 (D' 1 ) k - 1 = 2b (-a 1 ) k x, respectively Evaluating 
the averages of products of three multiple stochastic integrals and summing over all the 
triplets (n, m, I) contributing to the same order in r, as done for the second order coefficient, 
one gets D 3 = — c\{r) + c 3 (r)x 3 , as in Eq. §M$- 

Analogously, at fourth order, considering the relevant products of four multiple integrals, 
for 62 = 0, three types of contributions appear, yielding D^x, r) = do(r) + d2{r)x 2 + ^(r)^ 4 : , 
as in Eq. ([25]) . 

Let us recall that the averages of products of n multiple stochastic integrals appearing 
in the nt/i-order term of the coefficients can be expressed, in general, as multinomial terms, 
whose summation over all the products has the form nil k + fi2^ k + f^3^ k + ■■■■fi n n k (with 
rational fii) for the order k in r. Moreover, products of multiple stochastic integrals can 
be readily simplified, by means of useful relations between multiple Ito integrals [4]. Then, 
although at the cost of increasing the number of indices, the number of factors can be 
reduced. 
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